%% Michael D. Konig, Andrei Levchenko, Tim Rogers and Fabrizio Zilibotti
%% Aggregate fluctuations in adaptive production networks
%% Forthcoming in Proceedings of the National Academy of Sciences, 2022
 
%% Replication file (Matlab) for Figure 3, main text.

%% Requirements: 

% "rewiring_times_us_manuf.csv" -> CSV file with the following colums (US manufacturing firms only): replacement_year; lost_year; time; censored
% "rewiring_times_us.csv" -> CSV file with the following colums (US firms only): replacement_year; lost_year; time; censored
% "rewiring_times_manuf.csv" -> CSV file with the following colums (manufacturing firms only): replacement_year; lost_year; time; censored
% "rewiring_times.csv" -> CSV file with the following colums (full sample): replacement_year; lost_year; time; censored
% "rewiringtime_log_likelihood.m" -> Log-likelihood function for rewiring times.

clear all
close all

%% Parameters.
drop_non_us = true;
drop_non_manuf = false;
lb = []; %[0 0]; % Lower bound for bounded optimization algorithm.
ub = []; %[1 1]; % Upper bound for bounded optimization algorithm.

opt_alg = 'fminunc'; % choose from: 'nm' 'fmincon' 'fminunc' 'sa'
global options_fmincon;
global options_QP;
global options_NM;
global options_fminunc;
global options_SA;
options_fmincon = optimset('MaxIter', 1e10,'TolX',1e-5,'MaxFunEvals',10^10,'Display','off');
options_QP = optimset('Display','off');
options_NM = optimset('MaxIter', 1e10, 'TolX', 1e-5,'Display','off');
options_SA = optimset('MaxIter', 1e10, 'TolX', 1e-5,'Display','off');
options_fminunc = optimset('MaxIter', 1e10, 'TolX', 1e-5,'Display','off');

%% Load replacement-time data.
disp(['Load replacement-time data from file...'])
if(drop_non_us)
    if(drop_non_manuf)
        fid = fopen(['Data/rewiring_times_us_manuf.csv']); % replacement_year; lost_year; time; censored
    else
        fid = fopen(['Data/rewiring_times_us.csv']); % replacement_year; lost_year; time; censored
    end
else
    if(drop_non_manuf)
        fid = fopen(['Data/rewiring_times_manuf.csv']); % replacement_year; lost_year; time; censored
    else
        fid = fopen(['Data/rewiring_times.csv']); % replacement_year; lost_year; time; censored
    end
end
dat = textscan(fid,repmat('%s ',1,4), 'delimiter', ';', 'headerlines', 1);
fclose(fid);
rewiringtime = str2double(strtrim(dat{3}));
rewiringtime_censored = str2double(strtrim(dat{4}));
clear('dat')

W = nansum(rewiringtime);
U = nansum(1-rewiringtime_censored);
disp(['Theor. parameter estimate for gamma: ' num2str(U/W)])
disp(['Theor. estimated standard error for gamma: ' num2str(U/(W^2))])
theta_0 = U/W;

%% Compute maximum likelihood life-time estimates.
disp(['Compute maximum likelihood life-time estimates...'])
switch opt_alg
    
    case 'fmincon'
        
        disp(['using a constrained optimization algorithm (fmincon)...'])
        [theta_opt,fval,exitflag,output,lambda,grad,hessian] = fmincon(@(theta) -rewiringtime_log_likelihood(theta,rewiringtime,rewiringtime_censored),theta_0,[],[],[],[],lb,ub,[]);%,options_fmincon);
        
    case 'fminunc'
        
        disp(['using an unconstrained optimization algorithm (fminunc)...'])
        [theta_opt,fval,exitflag,output,grad,hessian] = fminunc(@(theta) -rewiringtime_log_likelihood(theta,rewiringtime,rewiringtime_censored),theta_0);%,options_fminunc);
        
    case 'nm'
        
        disp(['using an unconstrained Nelder-Mead algorithm...'])
        [theta_opt,fval] = fminsearch(@(theta) -rewiringtime_log_likelihood(theta,rewiringtime,rewiringtime_censored),theta_0);%,options_NM);
        
    case 'sa'
        
        disp(['using a constrained simulated annealing algorithm...'])
        [theta_opt,fval,exitFlag,out] = simulannealbnd(@(theta) -rewiringtime_log_likelihood(theta,rewiringtime,rewiringtime_censored),theta_0,lb,ub);%,options_SA);
        
    otherwise
        
        error('No optimization algorithm specified.')
end

disp(['Parameter estimate for gamma: ' num2str(theta_opt)])

if(strcmp(opt_alg,'fmincon') || strcmp(opt_alg,'fminunc'))
    err = sqrt(diag(inv(hessian)));
    disp(['Estimated standard error for gamma: ' num2str(err')])
end

z = theta_opt./err';
pvalue = 2*(1 - normcdf(z));
disp(['P-values for gamma: ' num2str(pvalue)])
